by R. Grothmann
Euler has a two routines to solve integer linear programs: intsimplex() and lpsolve(). intsimplex() is implemented in the Euler language, and uses the branch and bound method. lpsolve() loads the LPSOLVE library, which has been ported to Euler by Peter Notebeart.
This notebook demonstrates test cases.
Setup a test problem.
>A:=[10,8;9,11]; A:=1/A; ... b:=ones(rows(A),1); ... x:=0:0.1:15; y:=(b-A[:,1]*x)/A[:,2];
Plot the feasible restrictions.
>plot2d(x,y,a=0,b=12,c=0,d=12);
The function feasibleArea() computes the corners of the feasible set. This can be used to plot the area.
>xa:=feasibleArea(A,b); ... plot2d(xa[1],xa[2],filled=1,style="/",outline=1,add=1):
And this is the target function.
>c:=[1,1]
[1, 1]
We need <=.
>e:=-ones(size(b))
-1 -1
Solve this with the simplex method in Euler.
>{xs1,r,i}:=simplex(A,b,c,e,max=1); xs1,
7.10526 2.31579
The solution in fractional format.
>fraction xs1;
135/19 44/19
With intsimplex(), we find the optimal value for integer variables.
>intsimplex(A,b,c,max=1)
9 0
Or, if we require only y to be integer.
>fraction intsimplex(A,b,c,i=[2],max=1)
81/11 2
Let us add the solutions to the plot.
>plot2d(xs1[1],xs1[2],>points,>add); ... plot2d(9,0,>points,>add); ... plot2d(81/11,2,>points,>add):
Here is practical problem.
We take a square grid and want to select as many points as possible, with the restriction that at each point at most one of the neighbors (including the point itself) is chosen.
Euler has an incidence function, which can generate the incidence matrix for a graph, or the incidence matrix of a rectangle.
>load incidence;
We set up the incidence matrix, and add that the points are connected to themselves. The incidence matrix contains n squared entries such that
>n:=6; M:=id(n*n)+makeRectangleIncidence(n,n);
Add the restrictions that all points may be chosen only once.
>A:=M_id(n*n);
So we have the conditions
>b:=ones(rows(A),1);
The target function is the sum of x.
>c:=ones(1,cols(A));
All inequalities are <=.
>e:=-ones(rows(A),1);
The Euler Branch and Bound.
>xs1:=intsimplex(A,b,c,>max);
Here is the layout of the chosen nodes.
>format(4,0); (redim(xs1',[n,n])), longformat;
1 0 0 1 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 1 0 0 1 0 0 0
If you look ath the solution each entry in the grid does indeed have only one selected neighbor.
Note that we have n^2=36 variables and 72 restrictions. We can solve only a 6x6 grid fast enough with our elementary intsimplex() method.
Using LPSOLVE, it is possible to solve much larger problems.
>n:=12; M:=id(n*n)+makeRectangleIncidence(n,n); ... A:=M_id(n*n); ... b:=ones(rows(A),1); ... c:=ones(1,cols(A)); ... e:=-ones(rows(A),1); ... xs1:=ilpsolve(A,b,c,e,max=1); ... format(4,0); (redim(xs1',[n,n])), defformat;
1 0 0 1 0 0 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 1 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 1 0 0 0 0 1 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 1 0 0 0 1 0 0 1 0 0 0 1
Here is a random problem of the same kind.
We generate 100 points at random.
>n=100; x=random(n,1); y=random(n,1);
Then we generate a random incidence matrix with density 1/100.
>d=1/n; M=random(n,n)<d; M=max(M,M'); M=setdiag(M,0,1);
Compute the connected points and plot them.
>{xc,yc}=getGraph(x|y,M); ... plot2d(xc,yc,color=lightgray,<grid); ... plot2d(x,y,a=0,b=1,c=0,d=1,>points,>add):
Now we do our combinatorial optimization.
>A:=M_id(n); ... b:=ones(rows(A),1); ... c:=ones(1,cols(A)); ... e:=-ones(rows(A),1); ... xs1:=ilpsolve(A,b,c,e,max=1);
And plot those points in red, which are chosen by the algorithm.
>i=nonzeros(xs1'); ... plot2d(x[i],y[i],a=0,b=1,c=0,d=1,>points,color=red,style="[]#",>add):
The Knapsack problem is a famous combinatorial problem, where you have items with values and restrictions (usually the weight of the items), and a restricted space. You want to select items with as much value as possible.
The following example is from the Rosetta page. You have the following items.
>items=["map","compass","water","sandwich","glucose", ... "tin","banana","apple","cheese","beer","suntan cream", ... "camera","t-shirt","trousers","umbrella","waterproof trousers", ... "waterproof overclothes","note-case","sunglasses", ... "towel","socks","book"];
These items have the following weights each.
>ws = [9,13,153,50,15,68,27,39,23,52,11, ... 32,24,48,73,42,43,22,7,18,4,30];
And for your trip, they have the following values.
>vs = [150,35,200,160,60,45,60,40,30,10,70, ... 30,15,10,40,70,75,80,20,12,50,10];
You are not allowed to take more weight than a total of 400 with you, and you can only take one item of each kind.
The model is
This can immediately be translated into an integer linear program.
>A=ws_id(cols(ws)); b=[400]_ones(cols(vs),1); c=vs;
The matrix A contains the restriction in the first row, followed by rows of the type
>sol = intsimplex(A,b,c,eq=-1,>max,>check);
The solution is to take the following items.
>items[nonzeros(sol)]
map compass water sandwich glucose banana suntan cream waterproof trousers waterproof overclothes note-case sunglasses socks